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I. INTRODUCTION 

In recent experiments and numerical simulations it has been found that significant trans¬ 
port might be mediated by coherent structures such as streamers, blobs and vortices through 
the formation of rare avalanche-like events of large amplitude^ -8 . These events cause the 
probability distribution function (PDF) to deviate from a Gaussian profile on which the 
traditional mean field theory such as transport coefficients is based. Specifically, the PDF 
tails manifest the intermittent character of transport due to rare events of large amplitude 
that are often found to substantially differ from Gaussian distribution, although PDF cen¬ 
tres tend to be Gaussian. Therefore, a comprehensive predictive theory is called for in order 
to understand and subsequently improve intermittent transport features e.g. confinement 
degradation in tokamaks. 

Drift wave turbulence is known to generate zonal flows, which in turn inhibits the growth 
of turbulence and transport 1 ^. As such, zonal flows play an important role in fusion 

plasmas^ - —. 

In geophysical fluid dynamics, zonal flows are believed to cause a similar reduction in 
transport in atmospheres^ 1 ^ under certain limiting conditions^. This comes as no surprise, 
given the analogy between drift waves in the dissipationless limit and Rossby waves in 
nearly incompressible, shallow rotating fluids; both systems are described by the Charney- 
Hasegawa-Mima (CHM) equation—. 

Numerical studies where sheared flow is externally prescribed-^ - — lead to energetics 
that are qualitatively different from those obtained in the drift-wave/zonal-flow feedback 
mechanism^-. In the CHM model- 22- — sheared flow may be imposed by prescribing the 
background density gradient^. The known solutions of the CHM equation, namely bipolar 
and monopolar vortices^, form in the plasma. Here the term vortex is used to describe a 
localized extremum in the electrostatic potential that is evolved through the CHM equation. 
In this paper fluid simulations using the CHM equation are designed so that zonal flows 

are prescribed externally, 

moving in the poloidal direction and with various configurations in the radial direction. At 
initialisation bipolar vortices form, but only monopolar vortices survive due to the interaction 
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between vortices as well as the destructive effects of the sheared flow^. 


The fluid simulations produce quasi-stationary time series (poloidally averaged and sam¬ 
pled at different radial points) of the electrostatic potential and corresponding vorticity that 
describe the CHM flows^- 5,27,28 . We apply a standard Box-Jenkins modeling for each time 
series. This mathematical procedure effectively removes deterministic autocorrelations from 
the time series, allowing for the statistical interpretation of the stochastic residual part. In 
this particular case it turns out that an ARIMA(3,1,0) model (autoregressive integrated 
moving average)^ accurately describes the stochastic process. 


The stochastic residual of the time series of potential and vorticity exhibit Gaussian 
statistics or distributions with elevated exponential tails. We utilize analytical results from 
nonperturbative stochastic theory, the so-called instanton method^ - — for computing PDFs 
in turbulence as a comparison to the numerical data. The analytically derived PDFs are 
rather insensitive to the details of the linear physics of the system^ and thus display salient 
features of the nonlinear interactions. The numerically generated time traces are analysed 
using the ARIMA model and fitted with the analytical models accordingly. We find in the 
regions with strong nonlinear characteristics an emergent universal scaling of the PDF tails 
of exponential form ~ exp ( — const |0|) as suggested by recent theoretical work in Ref. 


36], 


38l . and 


40l relevant for the direct cascade dynamics. However, in many cases for the 


CHM zonal flows we find Gaussian PDFs in similarity to what was seen from the theoretical 


39 


, whereas for some mid radial points the system exhibit sub-exponential 


model in Ref. 

PDFs ~ exp ( — const \(f>\ a ) with a < 1, where the dynamics is strongly influenced by the 


zonal flow resulting in strong intermittency. 


The paper is organized as follows: Section [TT] introduces the CHM model and Section 1III1 
the statistical method used in the analysis. The numerical results are presented in Section 
HVl after which they are discussed and the paper is concluded by a summary of the main 
results. 
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II. CHARNEY-HASEGAWA-MIMA MODEL 


The CHM equation is solved in the Cartesian plane perpendicular to a constant magnetic 
field B = BqZ, with x and y being the radial and poloidal coordinates of toroidal geometry 
respectively. The plane is periodic in the y direction and finite in the x direction. The 
background electrostatic field 0 evolves through time on the (x, y ) plane through the CHM 

equation 22- — 


L n d 


fa - pMb - p. 


Ps C 


fa,pMd=°> 


where 


c s dt^ r> dy v * B 

are the Poission brackets. The diamagnetic velocity v* = u*y is defined as 


-l 


, - J_<9no _ c 2 s 

eB n 0 dx il cl 71 ’ 


( 1 ) 


( 2 ) 


with c the speed of light, T e the electron temperature, e the electron charge, c s the sound 
speed, Vt ci the ion cyclotron frequency and Uq the time-independent nonlinear background 
density. The characteristic length is the thermal Larrnor radius p s = c s /Vt ci and the char¬ 
acteristic time is taken as p s / max|u*| = (c s max \L~ X \)^ 1 . Equation (pQ) has two global 
invariants: the generalised energy W and the generalised enstrophy U. 

A second-order modified Euler predictor-corrector time scheme is used to solve equation 
(JTj). The periodic y direction is treated spectrally while the x direction as well as the non¬ 
linearity of equation (pQ) are finite differenced?^. The CHM equation’s nonlinear term is 
calculated using a conservative scheme for vector nonlinearities^. The x boundary condi¬ 
tions are fixed with 0 = 0 at all time. 

The numerical runs are initialized with a perturbation along the y direction consisting of 
many wavelengths, while the nonlinear density gradient L~ l = L“ 1 (a:) is prescribed and kept 
constant. The CHM equation produces a solution containing many pairs of bipolar vortices 
that evolve into larger monopolar vortices, the latter existing for most of the numerical run. 
In order to do a statistical analysis on these fluctuations, time series of the algebraic averages 
in the poloidal (y) direction of the normalised electrostatic potential ecf) / T e and normalised 
vorticity u;/f2 c ; are obtained at positions along the radial (x) axis, denoted as <fi(x,t) and 
u)(x,t). The normalised vorticity is obtained using 


co 

f - ri. 


= P 2 sV 2 


e<j> 


( 3 ) 
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III. STATISTICAL ANALYSIS 


In this paragraph we will quantify the intermittency in the simulated time series by 
computing the PDFs of the residuals or the stochastic component of the time traces and 
compare these with analytical predictions. Here, we briefly outline the implementation of 
the instanton method. For more details, the reader is referred to the existing literature^. 
In the instanton method the PDF tail is first formally expressed in terms of a path integral 
by utilizing the Gaussian statistics of the forcing, in a similar spirit as in Refs. |3l|, [34- 
371 . Here and throughout this paper, the term forcing is meant to describe the inherent 
unpredictability of the dynamics and will be_assumed to be Gaussian for simplicity. A 

The integral in the action ( S \) in the path 


general class of solutions is presented in Ref. 


36 


integral is evaluated using the saddle-point method in the limit A —> oo representing the tail 
values. The parameter A is proportional to some power of the quantity of interest such as the 
potential or flux. In mathematical terms, this corresponds to evaluating the integral along 
an optimum path among all possible paths or functional values. The instanton is localized in 
time, existing during the formation of coherent structure. The saddle-point solution of the 
dynamical variable 0(x, t) of the form 0(x, t) = F{t) r ift(x.) is called an instanton if F(t) = 0 
at t = —oo and F(t) ^ 0 at t = 0. Note that, the function ift( x ) here represents the 
spatial form of the coherent structure. Thus, the intermittent character of the transport 
consisting of bursty events can be described by the creation of the coherent structures. The 
dynamical system with a stochastic forcing is enforced to be satisfied by introducing a larger 
state space involving a conjugate variable (ft*, whereby (ft and (ft* constitute an uncertainty 
relation. Furthermore, (ft* acts as a mediator between the observables (potential or vorticity) 
and instantons (physical variables) through stochastic forcing. Based on the assumption that 
the total PDF can be characterized by an exponential form and that it is symmetric around 
the mean value fi, the expression 


ex P { —-^10 — Ml X }, 


(4) 


is found, where the potential (ft plays the role of the stochastic variable, with P((ft) deter¬ 
mining its statistical properties. Here b is a constant containing the physical properties of 
the system. Using the instanton method we find different statistics in different situations. 
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In a vorticity conserving system the intermittent properties of the time series in simulations 
are attributed to rare events of modon like structures that have a simplified response for the 
vorticity, 

V 2 ± (f) = — k±<j> + qx. (5) 

Here rj = 1+(1 —k 2 )U and the vortex speed is U. In this situation it has been predicted^ 1 ^ 
that the system has exponential tails in the direct cascade, exp ( — const | y|) ~ exp ( 
exp ( — const \k\(f)\), indicating a value of x = 1.0 as in Ref. 


const |V/ 
the References 


38 and 


40. In 


M 


35 


, and |37| the statistics of the momentum flux is found to be a stretched 


exponential with x = 3/2. However, when the nonlinear interactions are weak, as well as in 


the case o: 

on in Ref. 


an imposed zonal flow we find Gaussian statistics where x = 2 as is elucidated 


39. In the analysis we will make use of different types of distributions to retro-fit 


the PDFs of simulation results mainly using the Laplace distribution (x = 1-0) and the 
Gaussian distribution (x = 2.0). 

We focus on the time traces (averaged in the ^/-direction) at five equidistant radial points 
located at x — —18.9,-9.5,0.0,9.5,18.9 (in units of p s ). Each set of data describes the 
time evolution of the potential and vorticity to which we apply a standard Box-Jenkins 
modeling^. This mathematical procedure effectively removes deterministic autocorrelations 
from the system, allowing for the statistical interpretation of the residual part, which a 
posteriori turns out to be relevant for comparison with the analytical theory. In our set-up, 
it turns out that an ARIMA(3,1,0) model accurately describes the stochastic procedure, in 
that, one can express the (differenced) potential time trace in the form 


4>t +1 — «1 (pt + 4>t -1 + &3 <Pt-2 + (presit) 


( 6 ) 


where the fitted coefficients 01 , 02,03 describe the deterministic component and (p res is the 
residual part (noise or stochastic component). I 11 the time traces the mean values differ 
by several orders of magnitude and a convenient way of starting the comparison between 
different cases is to apply rescaling of the data. I 11 the rescaling procedure we multiply the 
original time trace with a constant factor, thus the mean and variance values are directly 
affected. However, the skewness and kurtosis are kept constant by construction. The benefits 
gained from rescaling are that we may compare a large number of different cases at different 
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radial points and that the tails are retained and the ARIMA model is preserved, thus in this 
sense the original and rescaled data is statistically equivalent. The original simulation data 
sets are down-sampled and consists of typically 5 x 10 4 entries. 

IV. RESULTS 

In this section the numerical results from all the different stream configurations are pre¬ 
sented in tandem with the statistical analysis. 

Throughout the paper the simulation plane has dimensions — L x < x < L x and 0 < y < 
L y where L x = 23.5 p s and L y = 23Ap s . The characteristic length is p s = 0.42 cm. 

A. One stream 

The constant background density gradient (Figure [l]) generates one stream centred at 
x = 0 and with velocity max |r>*| in the negative poloidal (y) direction. The flow is zero for 
|x| > 15/v The characteristic length and time in Figure [Ogive nrax |r;*| = 10 5 cm s _1 . 

At initialisation many vortices form that coalesce into large monopolar vortices, their 
widths determined by the width of the stream and with alternating polarities (Figure [2]). 
These monopolar vortices are dragged by the diamagnetic flow and move at a speed of 
0.2 max |n*| in the negative y direction, which is the direction of the diamagnetic velocity v* 
(Figured]). Smaller vortices that earlier in the evolution moved outside the stream, move 
very slowly in random directions under the influence of the large vortices inside the stream, 
as can be seen at positions x = ±20p s . 

After initialisation the generalised energy and enstrophy change significantly as vortices 
merge but once the large vortices have formed, from time 0.5 x 10 4 (c s max |-U” 1 1) —1 onwards, 
these quantities are relatively stable (Figure [3] and Table [I]). Small changes in energy con¬ 
servation are reflected in small changes in potential amplitudes. The change in enstrophy 
conservation is mirrored in changes in the vorticity of the fluctuations. 

In the statistical analysis, the higher moments of the distribution function may reveal 
important features of the statistics and the underlying dynamics of the process. Here, in the 
statistical analysis we will in addition consider the kurtosis and skewness. We define the kur- 
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FIG. 1. Constant x profile of the normalised diamagnetic velocity u*. The characteristic time is 
(< c s L - 1 )" 1 = 4.2 x 10~ 6 s. 


tosis as the fourth moment divided by the square of the second moment kurtosis = 7714 /m|, 
note that sometimes 3 is subtracted from the kurtosis yielding a zero kurtosis for a Gaus¬ 
sian distribution. A high value of the kurtosis is a key mark for a heavy tailed distribution 
which is flat at the centre. The skewness is defined as the third moment normalized by the 
3/2-power of the second moment, skewness = mz/rrn/ 2 and describes the asymmetry of the 
PDF around its mean value. 



FIG. 2. The normalised potential e4>/T e and normalised vorticity w/D c j at 23921.47 characteristic 
time units for the one stream. The minimum and maximum values of e<//T e are ±0.019 and the 
extrema for cu/f \ c i are -0.105 and 0.136. Maximum is white and minimum black in both plots. The 
vortices move in the negative y direction at a speed of 0.2 max |v*|. 
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conservation 



FIG. 3. Conservation of the generalised energy (dashed line) and enstrophy (dotted line) for the 
one stream. The change in generalised energy is AE/Eq = (Eq — E\)/Eq up to the end of the run 
and the change in generalised enstrophy is AU /Uq = (Uq — U\)/Uq. The notation E 0 and Uq is 
used for the initial values and E\ and U i for the final values. Table |U shows the conservation from 
times 0.48 x 10 4 and 0.96 x 10 4 (c s max |Zy~ 4 1) —1 onwards. 

A E/Et = (Et - E\)/E t AU/Ut = (U t - lh)/U t 

t = 0.2 t = 0.4 t = 0.2 t = 0.4 

One stream 0.026 0.171 —0.048 —0.068 


TABLE I. Conservation of the generalised energy E and generalised enstrophy U for the one stream. 
The time parameter t is scaled such that t = 1 is the end of the simulation. 

In Figure [1 the kurtosis along the x direction for the original time series is compared 
to the ARIMA modeled residual stochastic part of the time series. At some negative x 
locations distributions with elevated tails are found in the potential for the original time 
traces however this behaviour is not found for the vorticity. Furthermore, comparing the 
kurtosis of the potential and vorticity of the ARIMA modeled time traces the region with 
elevated tails coincide, as an indication of Eq. (J5]). 

Figure [5] displays the PDFs at several positions marked by black lines in Figure [4] along 
the x direction. The time traces are normalized in order to be able to show several PDFs 
in the same graph according to PDF(cj)) = PDF[(<p res — /x)//x], where /i is the mean value 


9 









Statistical properties of CHM zonal flows 


and <fi res is the residual after the ARIMA process. An analogous definition is adopted 
for the vorticity a). Here it can be seen that the PDFs of potential and kurtosis have 
exponential tails at the positions with higher kurtosis compared to the middle region which 
seems to exhibit Gaussian statistics. Note at some radial positions we find very large values of 
kurtosis signifying distributions with heavy tails and y < 1. Moreover, the PDFs are nearly 
symmetric yielding small values of the skewness measure. In particular, Figure [2] shows that 
the time evolution of the electrostatic potential exhibits a structure of alternating negative 
and positive polarity vortices along the poloidal (or y ) direction, indicating that the poloidal 
average for the nonlinear term is weak. At the same time the vortices are not symmetric 
about the x = 0 axis and these asymmetries result in the vorticity-likc statistics. 




FIG. 4. The kurtosis of the potential and vorticity time traces along the x direction in the one 
stream case are shown for the original time traces (graph on the left) compared to the ARIMA 
modeled time traces (graph on the right). 
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FIG. 5. Numerical PDFs along the x direction of the potential (graph on the left) and vorticity 
(graph on the right) time traces in the one stream case. The various x positions sampled are 
indicated by black vertical lines in Figured 
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B. Two adjacent anti-parallel streams 


The background density gradient (Figure [6]) generates two streams flowing anti-parallel 
to each other with max |u*| centred at x = ±3.5p s . The two streams flow adjacent to each 
other so that the edges of the streams are at x = ±14 p s and x = 0, where u* becomes zero. 
Two cases with different max |v*| are studied: a slow-flowing and a fast-flowing system. For 
the slow-flowing streams the characteristic length and time give max|v*| = 11.2 x 10 4 cm 
s _1 . The fast-flowing streams have a max |r;*| of 3.6 times that of the slow-flowing streams. 
Consequently the diamagnetic velocity gradient is larger for the fast-flowing streams. 


As in the case of one stream (Section IIV AD . the initialisation creates many vortices that 
merge to form monopolar vortices moving in each of the streams. In the case of the slow- 
flowing streams, a smaller gradient in the diamagnetic velocity shear exists between the two 
streams, compared to the shear in the case of the fast-flowing streams. As a result, there 
are more interactions between the vortices in the anti-parallel streams in the case of the 
slow-flowing streams. This manifests itself in Figure 0 by the fact that the two vortices 
situated inside the different streams move in the same direction, while the two vortices in 
Figure 0 move in opposite directions. 


Positive vortices do not move in the direction of the stream where they are situated. 
Instead, they tend to move in the direction of the stream that is sampled by the edge of the 
vortex. This is clearly shown in Figures 0 and 0 The negative vortices are dragged in the 
direction of their neighbouring positive vortices, shown by the vortex at x = 0 in Figure 0 
and the one at x = — lip,, in Figure 0 When the negative vortices are far enough from a 
positive vortex, their movement is determined by the flow direction sampled by their edge, 
as shown by the vortex at x = 15 p s in Figure 0 When vortices are far removed from the 
central flow streams and the vortices residing there, they exhibit small random movements, 
shown by the positive vortex at x = —15 p s in Figure 0 and the negative vortex at x = 16 p s 
in Figure 0 


The generalised energy and enstrophy conservation follow the same pattern as for one 
stream (Section IIV AD : after initialisation and large changes in energy and enstrophy the 
simulation produces a solution with less changes in these quantities. This is shown in Table 
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FIG. 6. Constant x profile of diamagnetic velocity n* normalised to max|a*| for two anti-parallel 
streams. For slow flows the characteristic time is (c s L~ 1 )~ l = 3.8 x 10 - ' 5 s and for fast flows it is 
(c.L- 1 )- 1 = 10- 5 s. 


A E/Et = (E t 

- E^/Et 

A U/U t = ( U t 

- Ui)/U t 

t = 0.2 

t = 0.4 

t = 0.2 

t = 0.4 

Slow-flowing streams 0.50 

0.29 

0.63 

0.37 

Fast-flowing streams 0.41 

0.26 

0.52 

0.30 


TABLE II. Conservation of the generalised energy E and generalised enstrophy U for two anti¬ 
parallel streams. The time parameter t is scaled such that t = 1 is the end of each simulation. 

nn where the change measured from 40% of the duration to the end of the time-line is 
approximately half the change measured from 20% to the end. The energy and enstrophy 
change throughout the run as the two adjacent streams influence each other. There is less 
interaction between the two streams when the diamagnetic velocity gradient between them 
is larger. This manifests in slightly better conservation in the faster streams (Table |TT|) . 

As a comparison to the statistical analysis of the one stream case, we will consider the 
time evolution of the stochastic part of the electrostatic potential and vorticity for both slow- 
flowing and fast-flowing anti-parallel streams. In Figure (9]( a), the kurtosis of the potential 
and vorticity time traces are displayed corresponding to the simulation of the slow flow 
presented in Figure [7] 
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We find similarly good correspondence in kurtosis profiles between the potential and 
vorticity as was found in the case of one stream (Section 11V AD for the stochastic residual 
part whereas for the original time traces no such correspondence could be found. We have 
omitted the figures of the kurtosis of original time series due to space limitations here and 
in the rest of the paper, since these do not provide any additional useful information. 

In Figure ITOl and fill the numerical PDFs of the stochastic part <fi res are shown for the 
CHM simulations of Figured We find that at the edge of the stream the PDFs are close to 
Gaussian whereas at the stream centre other nonlinear features can be found. Here we will 
utilize Eq. ([4]) in the following cases: Laplacian distribution denotes the analytical model 
for x — 1-0, whereas a Gaussian PDF is represented by x — 2.0. The appearance of a 
Laplacian distribution at the stream (Figure ITTl) is suggestive of a vorticity conserving non¬ 
linear system—’-^ whereas a Gaussian distribution is likely for a weakly nonlinear system—, 
or dynamics impeded by a strong zonal flow^. 

e0/T a at time 1055.88 ( c s L n _1 ) _1 

20 

15 

10 

5 

0 

-20 -10 0 10 20 
X ( P* ) 

FIG. 7. The normalised potential e</>/T e and normalised vorticity u/Qd at 1055.88 characteristic 
time units for the slow-flowing streams. The minimum and maximum values of ecj)/T e are —1.7 x 
10~ 2 and 1.4 X 10 2 . The extrema for uj/£l c i are —9.3 x 10 2 and 0.13. Maximum is white and 
minimum black in both plots. The positive vortex at x = 5 p s moves in the negative y direction 
at an approximate speed 3.4 max |a*|, while the negative vortex at x = 0 moves at approximate 
2 max 1| in the negative y direction. It accelerates and decelerates every time the positive vortex 
passes it on its right hand side. The positive vortex at x = —15 p s drifts slowly and randomly, 
while the negative vortex at x = 15 p s moves at speed 0.4 max |a*| in the positive y direction. 



oj/0 c; at time 1055.88 ( c s L n ) 



-20 -10 o 10 20 

X ( A, ) 
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e0/T e at time 3847.67 ( c s L n 1 ) 1 




-20 -10 0 10 20 
X ( Ps ) 


FIG. 8. The normalised potential e</>/T e and normalised vorticity u/Sld at 3847.67 characteristic 
time units for the fast-flowing streams. The minimum and maximum values of e(j)/T e are — 1.8x 10 -2 
and 10~ 2 , while the extrema for c d/Q c i are —3.9 x 10” 2 and 0.18. Maximum is white and minimum 
black in both plots. The negative vortex at x = —11 p s moves at speed 0.9 max |u*| and the positive 
vortex at x = —4 p s moves at max|u*|, both in the positive y direction. The positive vortex at 
x = 6 p s moves at speed 0.4 max |u*| in the negative y direction while the negative vortex at x = 16 p s 
undergoes small and slow random movements. 



FIG. 9. The kurtosis of the potential and vorticity of the ARIMA modeled time traces along the 
x direction are shown for (a) the slow-flowing and (b) the fast-flowing anti-parallel cases. 


In Figure |^b) the kurtosis of the time series generated by the fast-flowing anti-parallel 
streams corresponding to the simulation presented in Figure [8j The ARIMA modeled time 
trace of potential and vorticity follow each other reasonably well compared to the original 
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FIG. 10. The PDFs of the potential (graph on left) and vorticity (graph on right) time traces 
along the x direction in the slow-flowing anti-parallel case. The various x positions sampled are 
indicated by black lines in Figure [9|a). 

time traces. Also here the skewness is small and the PDFs are Gaussian or exponential. 


16 
























Statistical properties of CHM zonal flows 





FIG. If. Fits to the two salient types of PDFs, the Gaussian and the Laplacian (exponential) 
distributions. 
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C. Two streams with space between them 

Figure fl2l shows the background density gradient generating the two cases studied here: 
two streams flowing parallel to each other and two streams flowing anti-parallel to each 
other. For both cases there exists a space between the streams. The parallel streams are 
located at positions x G [—14, —4] and x G [4,14] with max |u*| at x = ±9 p s and v* = 0 
outside these intervals. The anti-parallel streams are located at positions x G [—16, —6] and 
x G [6,16] with max |v*| at x — ±llp s and u* = 0 outside these intervals. The characteristic 
length and time are the same in both cases, giving max |u*| = 4.7 x 10 4 cm s _1 . 

Initial pairs of bipolar vortices merge to form monopolar vortices. In the case of two 
parallel flowing streams, all the vortices move in the flow direction of the two streams (Figure 
[T3l) . These vortices exhibit the same behaviour as those in Section IIVBI namely they move 
in the direction of the flow sampled by their edges. Vortices outside the two streams are not 
influenced by the flow direction of the streams, e.g., the vortices at x = ±20 p s in Figure fl3l 

In the case of the two anti-parallel flowing streams, monopolar vortices form inside the 
streams from bipolar vortices after initialisation. These vortices move in the opposite di¬ 
rection of the stream in which they reside. Monopolar vortices forming between the two 
anti-parallel streams at positions x > 0 move in the flow direction their edges sample, i.e., 
the negative y direction. At some stage during the simulation these vortices migrate in the 
negative x direction, moving through the stream situated at x — —10p s and in the process 
destroying the vortices in that stream. The result is shown in Figure ITT where one vortex, 
situated at x — —p s , moves between the streams in the negative y direction while another 
moves in the opposite direction of the stream flow at x — 8 p s . 

Table IIIII shows that the generalised energy conservation is good for the two streams 
flowing in the same direction, but bad for the anti-parallel streams. This shows that the 
amplitude of the fluctuations decrease more in the latter case. There is no discernible pattern 
in the conservation of the generalised enstrophy. 

Figure fl5l presents the kurtosis generated with two imposed parallel or anti-parallel flows 
with a spacing between them. The statistics are obtained from the time evolution of the 
electrostatic potential in Figures [13] and [T4] for the parallel and anti-parallel flows respec- 
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Parallel Flows 



X ( Ps ) 


Anti —Parallel Flows 



FIG. 12. Constant x profiles of diamagnetic velocity u* normalised to max |u*. [ for two streams 
with a space between them. The characteristic time is (c^L” 1 ) -1 = 8.9 x ICC 6 s for both parallel 
and anti-parallel streams. 



FIG. 13. The normalised potential e<f>/T e and normalised vorticity ca/f id at 1124.14 characteristic 
time units for two parallel streams. The minimum and maximum values of e</>/T e are —2.3 x 10 -2 
and 2.4 x 10~ 2 , while the extrema for ca/f2 c j are 0.18 and 0.2. Maximum is white and minimum 
black in both plots. The negative vortex at x = —3 p s moves at speed 0.7 max |u*| and the positive 
vortex at x = 6 p s moves at 1.2max |u*|, both in the negative y direction. 

tivcly. Note that the kurtosis of time traces of potential and vorticity follow each other 
reasonably well compared to the original time traces. 
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FIG. 14. The normalised potential e0/T e and normalised vorticity oj/Cld at 2273.71 characteristic 
time units for anti-parallel streams. The minimum and maximum values of e0/T e are —1.8 x 10“ 2 
and 2.5 x 10 —2 , while the extrema for co/Q C i are -0.25 and 5.2 x 10~ 2 . Maximum is white and 
minimum black in both plots. The negative vortex at x = — p s moves at speed 1.4 max |u*| in the 
negative y direction and the positive vortex at x = 8 p s moves at speed 0.6 max |n*| in the positive 
y direction. The positive vortex at x = —21 p s hardly moves at all. 


A E/E t = (E t 

- E^/Et 

A U/U t = ( U t 

- Ui)/U t 

t = 0.2 

t = 0.4 

t = 0.2 

t = 0.4 

Parallel streams 0.03 

-0.003 

0.007 

-0.21 

Anti-parallel streams 0.32 

0.160 

0.274 

0.02 


TABLE III. Conservation of the generalised energy E and generalised enstrophy U for two streams 
with a space between them. The time parameter t is scaled such that t = 1 is the end of each 
simulation. 
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FIG. 15. The kurtosis of the potential and vorticity of the ARIMA modeled time traces along the 
x direction are shown in the (a) parallel and (b) anti-parallel cases. 
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V. DISCUSSION AND SUMMARY 

In this paper simulations of unforced Charney-Hasegawa-Mima (CHM) flows have been 
performed where we have imposed constant background density profiles on the CHM equa¬ 
tion. In this manner we generated a single zonal flow, anti-parallel flows with slow and fast 
moving plasma adjacent to each other, as well as parallel and anti-parallel flows with a space 
between the flow streams. 

In the single flow stream a chain of monopolar vortices form along the stream with their 
widths determined by the stream width and with polarities opposite from their neighbours. 
The vortex chain moves in the flow direction. When two anti-parallel zonal flows are placed 
adjacent to each other, monopolar vortices form in each stream but they are affected by 
the neighbouring stream and vortices in their close proximity. The smaller the flow velocity 
gradients are, the more interaction occurs. Vortices move in the direction their edges sample. 
In addition, negative polarity vortices in the presence of positive vortices are dragged along 
by the positive vortex. When two streams are placed a distance from each other, monopolar 
vortices form inside as well as between the streams. All vortices move in the flow direction 
when the streams flow in parallel. In the case of anti-parallel flow the vortex movement is 
in the opposite direction of the residing flow velocities. Vortex movement between the two 
flows are dictated by the flow direction the vortex edge samples. 

After initialisation bipolar vortices form that merge and destroy each other until only 
monopolar vortices survive, which leads to enstrophy changes during the initial phase. The 
amplitudes of the vortices also change to a stable value during this phase, hence an initial 
change in energy. Once the initial phase is over, the energy and enstrophy conservation is 
good for the one stream. In the case of the two stream scenarios the energy and enstrophy 
are conserved to a lesser extent in some numerical runs and to a larger extent in others. 
This is because vortices keep on being created and destroyed as they interact in a random 
manner with each other and the sheared flow. 

We have sampled time series at points in the radial direction (poloidally averaged) of the 
electrostatic potential and corresponding vorticity generated by the simulations. Our aim is 
to evaluate the intermittent characteristics of the time series by using a standard Box-Jenkins 
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modeling. This mathematical procedure effectively removes deterministic autocorrelations 
from the time series, allowing for the statistical interpretation of the stochastic residual part. 
The numerically generated time traces are compared with predictions from a nonperturbative 
theory, the so-called instanton method for computing probability density functions (PDFs) 
in turbulence. More specifically the numerically generated time traces are analysed using the 
ARIMA model and fitted with analytical models accordingly. In the simulations presented 
here we find that an ARIMA(3,1,0) model presents an adequate description of the stochastic 
process. 

The time series of the ARIMA modeled stochastic residual of cases described in section 
[TV](A,B,C) of the potential and vorticity exhibit in general a uni-modal PDF with Gaussian 
features or a PDF with exponential tails. In summary, the PDF ~ exp ( — const |0| x ) with 
X = 2.0 or exponential statistics are found with (y = 1.0). The different configurations are 
represented by an imposed slow and a fast zonal flow as well as parallel and anti-parallel 
flows. The rationale for using the ARIMA model is to uncover the stochastic process hidden 
in the numerically generated time trace. The ARIMA process is specifically designed to 
identify correlations in a time trace by utilizing a differencing procedure and thus provides 
us with an efficient model how to subtract the correlations in time from the full signal. The 
objective of this work is to identify the particular stochastic process that is generating the 
time trace. The one important restriction of this procedure is that the stochastic process has 
to be stationary with respect to the mean and variance. For an arbitrary stochastic process, 
there exist formal tests to check whether stationarity holds. In this particular example of 
CHM zonal flows with a constant background density profile generating the process, this is 
fulfilled except for a short interval at the start of the simulation. 

One non-trivial aspect of the statistics of the electrostatic potential in the CHM system is 
conservation of energy and enstrophy that may lead to a Gibbsian equilibrium distribution. 
However with the imposed shear flow through the varying density gradient, non-vanishing 
triad interactions redistribute energy between the modes, contributing to a situation where 
the PDFs deviate from the Gaussian form at some radial locations. For a general discussion 
see Kronnnes^k The cascade processes or redistribution energy and enstrophy can also be 
directly seen in the simulations as deviations in the conserved quantities. 
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To this end, in general we find an emergent universal exponential scaling of the distri¬ 
bution functions (Laplace distribution with % = 1.0) that accurately describes the statistics 
of the time series of the electric potential and vorticity. Analysing the profiles along the 
x coordinate of the kurtosis of the potential and vorticity time series we find striking sim¬ 
ilarity suggestive of the relation in Eq. ()5|) . Exceptions to the named distributions occur 
where strong nonlinear interactions are present in the dynamics where the PDFs are sub¬ 
exponential with y < 1.0 and high values of the normalized fourth moment (kurtosis) is 
found. 
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